We construct a statistical model of crystal melting to count BPS bound states of D0 and D2 branes on a single D6 brane wrapping an arbitrary toric Calabi-Yau threefold. The three-dimensional crystalline structure is determined by the quiver diagram and the brane tiling which characterize the low energy effective theory of D branes. The crystal is composed of atoms of different colors, each of which corresponds to a node of the quiver diagram, and the chemical bond is dictated by the arrows of the quiver diagram. BPS states are constructed by removing atoms from the crystal. This generalizes the earlier results on the BPS state counting to an arbitrary non-compact toric Calabi-Yau manifold. We point out that a proper understanding of the relation between the topological string theory and the crystal melting involves the wall crossing in the Donaldson-Thomas theory.